\defmodule {GammaProcessPCABridge}

Same as \class{GammaProcessPCA}, but the generated uniforms
correspond to a bridge transformation of the \class{BrownianMotionPCA}
instead of a sequential transformation.

\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        GammaProcessPCABridge
 * Description:
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       Jean-Sebastien Parent and Maxime Dion
 * @since        july 2008

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.stochprocess;\begin{hide}
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;

\end{hide}

public class GammaProcessPCABridge extends GammaProcessPCA \begin{hide} {
    protected BrownianMotionBridge BMBridge;
    protected double       mu2OverNu,
                           mu2dTOverNu;
    protected double[]     bMu2dtOverNuL,  // For precomputations for G Bridge
                           bMu2dtOverNuR;
    protected int[]        wIndexList;
\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public GammaProcessPCABridge (double s0, double mu, double nu,
                                 RandomStream stream)\begin{hide} {
        super (s0, mu, nu,  stream);
        this.BMBridge = new BrownianMotionBridge(0.0, 0.0, Math.sqrt(nu), stream);
    }\end{hide}
\end{code}
\begin{tabb} Constructs a new \texttt{GammaProcessPCABridge} with parameters
$\mu = \texttt{mu}$, $\nu = \texttt{nu}$ and initial value $S(t_{0}) = \texttt{s0}$.
The same stream is used to generate the gamma and beta random numbers.  All
these numbers are generated by inversion in the following order:
the first uniform random number generated
is used for the gamma and the other $d-1$ for the beta's.

\end{tabb}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code}\begin{hide}

    public double[] generatePath (double[] uniform01) {
	// uniformsV[] of size d+1, but element 0 never used.
        double[] uniformsV = new double[d+1];

	// generate BrownianMotion PCA path
        double[] BMPCApath = BMPCA.generatePath(uniform01);
        int oldIndexL;
        int newIndex;
        int oldIndexR;
        double temp, y;
    // Transform BMPCA points to uniforms using an inverse bridge.
        for (int j = 0; j < 3*(d-1); j+=3) {
            oldIndexL   = BMBridge.wIndexList[j];
            newIndex    = BMBridge.wIndexList[j + 1];
            oldIndexR   = BMBridge.wIndexList[j + 2];

            temp = BMPCApath[newIndex] - BMPCApath[oldIndexL];
            temp -= (BMPCApath[oldIndexR] - BMPCApath[oldIndexL]) * BMBridge.wMuDt[newIndex];
            temp /= BMBridge.wSqrtDt[newIndex];
            uniformsV[newIndex] = NormalDist.cdf01( temp );
        }
	double dT = BMPCA.t[d] - BMPCA.t[0];
	uniformsV[d] = NormalDist.cdf01( ( BMPCApath[d] - BMPCApath[0] - BMPCA.mu*dT )/
					 ( BMPCA.sigma * Math.sqrt(dT) ) );


	// Reconstruct the bridge for the GammaProcess from the Brownian uniforms.
        path[0] = x0;
        path[d] = x0 + GammaDist.inverseF(mu2dTOverNu, muOverNu, 10, uniformsV[d]);
        for (int j = 0; j < 3*(d-1); j+=3) {
            oldIndexL   = wIndexList[j];
            newIndex    = wIndexList[j + 1];
            oldIndexR   = wIndexList[j + 2];

            y =  BetaDist.inverseF(bMu2dtOverNuL[newIndex],  bMu2dtOverNuR[newIndex], 8, uniformsV[newIndex]);

            path[newIndex] = path[oldIndexL] +
		(path[oldIndexR] - path[oldIndexL]) * y;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


    public double[] generatePath()  {
        double[] u = new double[d];
        for(int i =0; i < d; i++)
            u[i] = stream.nextDouble();
        return generatePath(u);
    }


    public void setParams (double s0, double mu, double nu) {
        super.setParams(s0, mu, nu);
        BMBridge.setParams(0.0, 0.0, Math.sqrt(nu));
        BMPCA.setParams(0.0, 0.0, Math.sqrt(nu));
    }


    public void setObservationTimes (double[] t, int d) {
        super.setObservationTimes(t, d);
        BMBridge.setObservationTimes(t, d);
    }\end{hide}

    public BrownianMotionPCA getBMPCA() \begin{hide} {
        return BMPCA;
    }\end{hide}
\end{code}
\begin{tabb} Returns the inner \class{BrownianMotionPCA}.
\end{tabb}


\begin{code}
\begin{hide}
    protected void init() {
        super.init();
        if (observationTimesSet) {

        // Quantities for gamma bridge process
        bMu2dtOverNuL = new double[d+1];
        bMu2dtOverNuR = new double[d+1];
        wIndexList  = new int[3*d];

        int[] ptIndex = new int[d+1];
        int   indexCounter = 0;
        int   newIndex, oldLeft, oldRight;

        ptIndex[0] = 0;
        ptIndex[1] = d;

        mu2OverNu   = mu * mu / nu;
        mu2dTOverNu = mu2OverNu * (t[d] - t[0]);

        for (int powOfTwo = 1; powOfTwo <= d/2; powOfTwo *= 2) {
            /* Make room in the indexing array "ptIndex" */
            for (int j = powOfTwo; j >= 1; j--) { ptIndex[2*j] = ptIndex[j]; }

            /* Insert new indices and Calculate constants */
            for (int j = 1; j <= powOfTwo; j++) {
                oldLeft  = 2*j - 2;
                oldRight = 2*j;
                newIndex = (int) (0.5*(ptIndex[oldLeft] + ptIndex[oldRight]));

                bMu2dtOverNuL[newIndex] = mu * mu
                                   * (t[newIndex] - t[ptIndex[oldLeft]]) / nu;
                bMu2dtOverNuR[newIndex] = mu * mu
                                  * (t[ptIndex[oldRight]] - t[newIndex]) / nu;

                ptIndex[oldLeft + 1]       = newIndex;
                wIndexList[indexCounter]   = ptIndex[oldLeft];
                wIndexList[indexCounter+1] = newIndex;
                wIndexList[indexCounter+2] = ptIndex[oldRight];

                indexCounter += 3;
            }
        }
        /* Check if there are holes remaining and fill them */
        for (int k = 1; k < d; k++) {
            if (ptIndex[k-1] + 1 < ptIndex[k]) {
            // there is a hole between (k-1) and k.

                bMu2dtOverNuL[ptIndex[k-1]+1] = mu * mu
                                  * (t[ptIndex[k-1]+1] - t[ptIndex[k-1]]) / nu;
                bMu2dtOverNuR[ptIndex[k-1]+1] = mu * mu
                                  * (t[ptIndex[k]] - t[ptIndex[k-1]+1]) / nu;

                wIndexList[indexCounter]   = ptIndex[k]-2;
                wIndexList[indexCounter+1] = ptIndex[k]-1;
                wIndexList[indexCounter+2] = ptIndex[k];
                indexCounter += 3;
            }
        }
        }
    }

}
\end{hide}
\end{code}
